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ABSTRACT 

We seek to develop an essentially analytical model for the evolution of Fanaroff-Riley 
Class II radio galaxies as they age individually and as their numbers vary with cos- 
mological epoch. Such modeling is required in order to probe in more detail the im- 
pact of radio galaxies on the growth of structures in the universe, which appears 
likely to have been quite significant at z > 1. In this first paper of a series we com- 
pare three rather sophisticated analytical mod els for the evolution of linear size and 
lobe power of FR II radio galaxies, those of IKaiser. Dennett-Thorpe fc Alexander! 
l|1997|h iBlundell. Rawlings fc WillottJ l|1999ft and iManolakou fc Kirkl l)2002ft . We p~ 
form multi-dimensional Monte Carlo simulations in order to compare the predictions of 
each model for radio powers, sizes, redshifts and spectral indices with data. The obser- 
vational samples used here are the low frequency radio surveys, 3CRR, 6CE and 7CRS, 
which are flux-limited and complete. We search for and describe the best parameters 
for each model, after doing statistical tests on them. We find that no existing model 
can give acceptable fits to all the properties of the surveys consi dered, although the 
IKaiser et alJ (I1997T) model give s better overall results than do the Manolak ou fc Kirkl 
(2002) or Blundell et al.l ljl999) models for most of the tests we performed. We suggest 
ways in which these models may be improved. 

Key words: galaxies: active - galaxies: evolution - intergalactic medium - large-scale 
structure of universe - methods: statistical - radio continuum: galaxies. 



1 INTRODUCTION 

Radio galaxies (RGs) with extended lobes on opposite sides 
of their nuclei, or the classical double s ources, constitute 
a sign ificant population of active galaxies. iFanaroff fc Rilevl 
( 1974) classified these objects as Class II (FR II) sources. 
These are the edge-brightened population of radio sources 
(typically with bright hotspots) , and also the more powerful 
ones, with luminosities Pi7smhz > 10 25 W Hz" 1 sr 

Flux limited samples indicate that the comoving den- 
sities of RGs were higher during the quasar era (i.e., 
between redshifts ~ 1.5 and 3) as compared to the 
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x-ray observations of powerful active ga lactic nuclei revea l 
a similar trend for the quasar era (e.g., lUeda et aljl2003ft . 
The star and galaxy fo rmation rate was also considerably 
higher in the quasar era. iLillv et ail lll996l) inferred that the 
observed luminosity density (and hence the star formation 
rate) of the universe in the UV, optical and near-infrared 
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increases markedly with redshift over < z < 1 . Simi larly, 
from Hubble Deep Field studies IConnollv et al.l il997l^ and 
iMadau. Pozzetti fc Dickinson! jl99ct) found a sharp rise in 
the comoving luminosity density and global star formation 
rate with redshift, finding that it peaked at z ~ 1.5, and 
decrease d monotonically at highe r z ou t to z ~ 3 — 4. More 
recently. iBouwens fc Illingworthl |2006) found an apparent 
decrease in the rest-frame UV luminosity function and the 
cosmic star formation rate density from the peak redshift 
of z ~ 3 upto z ~ 6 — 10. Studies made with t he Spitzer 
Space Telescope fe.g.. lPerez-Gonzalez et al]l2005l) also indi- 
cate that the infrared luminosity function and the cosmic 
star formation rate increase with redshift until the quasar 
era. 

Submillimeter surveys claimed that the c omoving lu- 
mino s ity density has a pea k at z ~ 2 — 5 fB lain et alJ 
Il999t lArchibald et~ai]|20Qj) . This redshift range is some- 
what higher than what optical surveys (possibly affected by 
dust obscu ration) infer. At the sa me time, a more recent sub- 
mm study jRawlings et al .120041) indicates no compelling ev- 
idence for the far infrared luminosity of radio sources to rise 
with redshift. 

The above observations have prompted investiga- 
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tions of the effect of RGs on cosmological evolution and 
distribution of large scale structures in the universe. 
Preliminary work indicates that RGs can have substantial 
impacts on the formation, distribution and evolution 
of gal axies and large scale s tructu res of the universe 
(e.g.. IGopal-Krishna fc Wiital l200ll. hereafter GKW01; 
iKronberg et all l200lt IGopal-Krishna fc Wiital 120031 
hereafter GKW03; lGopal-I^ rMvna :j _WiiJ 1 a L &: OstermarJ 

Ol hereafter GKWO; 
2004 hereafter G KWB: 
Levine fc Gnedi J2005H . 



Gopal-Krishna, Wiita & Barai 
iRawlines fc Jarvisl |2004 



One important aspect of this process is the role played 
by the huge expanding RG lobes in triggering extensive 
star formation in a multi-phase intergalactic medium. This 
idea has been discussed by several authors in order to ex- 
plain the alignment between lar ge scale optical emission 
and radio source direction (e.g ., iBeeelman fc C ioffi 1989; 
iDe Youndfl989l) . IChokshil dl997T) proposed that RG expan- 
sion could trigger much star formation in host galaxies. 
GKW01 stressed that RGs could impact a large fraction 
of the filamentary structures in which galaxies form, thus 
hastening their birth. Similar c onclusions were drawn from 
different lines of a r gume nt by IKronberg et al.l (l200lf) and 
iFurlanetto fc Loebl feOOlT) . Recently. ISilkl (l2005h also argued 
that efficient ultraluminous starbursts can be triggered by 
AGN jets. 

A very significant fraction of the volume of the uni- 
verse in which star formation has occurred was impinged 
upon by the growing radio lobes during the quasar era 
(GKW01 and references therein). When these radio lobes 
propagating through the protogalactic medium envelop 
cooler clumps of gas embedded within the hotter gas 
which fills most of the volume, the initial bow shock 
compression triggers large-scale star formation, which is 
sustained by the persistent overpressure from the engulf- 
ing radio cocoon. This cocoon pres sure is likely to be 
well a bove the equipartition estimate jBlundell fc Rawlingsl 
2000). This scenario is supported by many computa- 
tions, analyt ical (e.g., [Reesj Il989t iDalvl Il990 |). hydrodv- 
nami c al (e.g.. | De Y punp|l989: Mcllema . Kurk fc Rotteerind 
120021: IFragile et alJ 120041: ISaxton et alJ 120051) . and mag- 
netohydrodynamical (e.g., IFragile et al.l 120051) . This trig- 
gered star formation provides an explanation for much 
of the remarkable radio-optical al ignment effect exhib- 
ited by high-z radio galaxies (e.g.. [ McCarthy et alJ Il987t 
IChambers. Milev fc van Breugell Il988h . Additional sudd ort 
for jet or lobe-induced star formation comes from the 
Hub ble Space Telescope images of z ~ 1 radio galax- 
ies feest. Longair fc Rottgerind Il996f) and of so me ra- 
dio sources a t higher z (e.g. iBicknell et aljEoOoT) . Keck 



observations \~Dej et alj 
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1997) and sub-mm observations 



2005) of high z RGs also give ev- 



idence for this phenomenon. Clustered Lyman a emitters 
have been found at high redshifts (z ~ 2 — 5) close to RGs 
llVenemans et all2005l : Overzier et all20 05). indicating that 
RGs form in high density regions and could have significant 
impact by accelerating star formation. Deep optical HST 
imaging gives evidence of star formation and a starburst 
driv en superwind ind uced by AGN jet activity in a z = 4.1 
RG JZirm et al.ll2005l) . 

The expanding RG lobes also could have infused mag- 
netic fields of significant strengths (~ 10~ 8 Gauss, e.g., 



iRvu. Kang fc Biermannl Il998|) into the cosmic web por- 
tion of the IGM (GKW01: IKronberg et alJ l200ll : GKWO; 
GKWB). Evidence of substantial metalicity in underd ense 
regions of the IGM at z ~ 4 (e.g., ISchave et al"ll2003h re- 
quires a strong mechanism of spreading metals widely (met- 
alization) at early cosmic epochs. The huge radio lobes could 
contribute substantially to spreading metals into the IGM, 
by sweeping out the metal-rich ISM of some young galaxies 
which they encounter while expanding (GKW03, GKWB). 

Ascertaining the importance of these processes of star 
formation, magnetization and metalization via RGs requires 
addressing the question of what fraction of the relevant 
volume of the universe did the radio lobes occupy during 
the quasar era (GKW01). The "relevant universe" refers to 
the volume containing most of the baryons, the majority 
of which exist as a filamentary structure of warm/hot gas, 
the WHIM ( warm/hot intergalactic medium) with 10 5 < T 
< 10 7 K (e.g.. lCen fc Ostrikerll999l.l2006UDave et all200lTl . 
For the radio lobes to have an important role in impacting 
star formation and spreading magnetic fields and metals, 
they need to occupy a significant portion of this relevant 
volume of baryons, which, however, was a small fraction of 
the total volume of the universe during most of the quasar 
epoch. A prerequisite for a more accurate computation of 
this RG impacted volume is a good model of the evolution 
of radio sources, both individually and as a function of z. 

Many analytical models have been published which 
characterize radio so urces in terms of t h eir d ynamics 
and power evolution. iKaiser fc Alexander! il997l here- 
after KA) showed that the cocoons can have a self-similar 
growth. Although numer ical hydrodynamical studies (e.g., 
ICarvalho fc O'Deal I2OO2T1 indicate that radio source sizes 
grow in a more complex way than self-similar predictions, 
they are still reasonable approximations overall. The power 
evolutions in these models are dominated by adiabatic losses 
as the lobe expands, synchrotron radiation losses in the lobe 
magnetic field and inverse compton (IC) losses off the cosmic 
microwave background (CMB) photons. 

The three models of radio lobe power evolution with 
time wh ich are considered in detail in this paper are those 
given bv lKaise r 1 Dennett-Thorpe fc Alexander) <[l997l . here- 
after KDA) iBlundell. Rawlings fc Willottl (Il999l . hereafter 
BRW) and iManolakou fc Kirkl fcOOl hereafter MK). The 
source linear size evolution in BRW and MK essentially fol- 
low the KDA prescription. They differ in the way the rela- 
tivistic particles are injected from the jet to the lobe, and in 
treatments of loss terms and particle transport. So there are 
some significant differences in their predictions for observed 
powers (P) as functions of source size (D) and redshift (z). 

The simplest method to study the power evolution of 
RGs is to examine their radio power - linear size, or P-D, 
diagram. Thes e P-D tracks have b een used (KDA; MK; 
Machalsk i. Chvzv fc Jamrozvl l2004 a b) to look for consis- 
tency between data and models. These papers compare 
model tracks with P-D diagrams of observed radio sources 
to evaluate the qualitative success of the models. 

The innovative radio sky simulation prescription in 
BRW adds new dimensions to the observed parameter space. 
Using the R G redshift dist ri butio n estimated by BRW from 
the work of IWillott et alJ (1200 on the radio luminosity 
function (RLF) and any lobe power evolution model, one 
can get P, D, z, and spectral index a {P v <x v~ a ) values for 
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simulated model radio sources. The distributions of these 
simulated RGs can then be compared to observational data 
to test the success of the model. In BRW, slices through the 
[P, D, z, a]-space generated by their model are qualitatively 
compared with observations for two data sets (3CRR and 
7CRS); those authors claim good results, except for plots 
involving a. 

However, to properly claim success for a theoretical 
model a quantitative statistical test is required; we present 
some in this paper. A quantitative comparison of cosmolog- 
ical radio source evolution model predicti ons with an obser - 
vational data sa mple (the 3C data fro mjLain g et al.lll983l) 
has been done bv lKaiser fc Alexander! dl993) . Thev consid- 
ered a progenitor FR II source population being born over 
cosmic epochs, and evolving according to assumed distri- 
bution functions of the model parameters of the KDA and 
KA models. Constructing simulated samples, they then com- 
pared the models' predictions with observations. They used 
X 2 statistics in the [P — D] and [P— z] planes to constrain the 
models. However the binning they used was somewhat arbi- 
trary and the bins appear to be based on the concentration 
of sources in the observed [P — D — z] planes. 

Our approach (based on 1- and 2-dimensional 
Kolmogorov-Smirnov (KS) statistics and correlation coef- 
ficients) may be as good as can be done since we are deal- 
ing with source characteristics in four dimensions (P, D, 
z, a) and over three observational surveys (3CRR, 6CE and 
7CRS) with only a few hundred sources in total. We tried to 
perform multi-dimensional KS like tests (discussed in §5.2.1) 
but the limited sizes of the observational samples precluded 
any useful results from being obtained. 

In §2, we summarize the BRW simulation prescription. 
We apply this prescription to the KDA, BRW and MK mod- 
els in §3. In §4 we discuss the observational samples to which 
we will compare the model distributions, and describe how 
our multi-dimensional Monte Carlo simulations are done. 
We perform statistical tests comparing the distributions of 
radio source parameters predicted by each model and those 
of observational samples in §5. We vary the parameters of 
the models, aiming to find the parameters which give the 
best statistical fit for each model to all three surveys simul- 
taneously. A discussion and conclusions are given in §6 and 
§7, respectively. 



2 INITIAL POPULATION GENERATION 

We follow the prescription given in detail in BRW to gener- 
ate the initial radio source population. Here we summarize 
the initial distributions of source ages, redshifts and beam 
powers; these produce the redshift, beam power and the age 
at which each model RG will be intercepted by our light 
cone. This summary and update of the BRW prescription is 
necessary to define the model parameters. One key difference 
from BRW is that we assume a consensus cosmology, i.e., a 
flat univer se with Hp = 71 km s _1 Mpc -1 , fi m = 0.3 and 
JIa = 0.7 iSoergel et alJ 2003). The cosmolog ical equations 
are taken from Ca rroll. Press fc Turned \l992) and lPeacockl 
(1999). 

From some initial high redshift, z s tart, well above the 
peak of the RLF, sources are assumed to be born at an 
interval, AT 3 t ar t, which is short compared to their lifetimes. 



From the cosmology assumed the reds hifts are tra nslated to 
cosmic times (epochs) and vice versa jWeinberdligsgh . We 
use Zstart = 10 and take AT,t ort = 10 6 years, but the results 
should be insensitive to values of z 3ta rt > 6 and AT 3ta rt < 
10 7 years. 

After a source is born at a redshift Zbirth, its active 
lifetime is denoted as TuaxAge.- A default value of TuaxAge = 
5 x 10 s years is taken. This value is used by BRW, and 
more recent investi gations involving X-ray activity in AGN 
( Barge r et al Il200ll) . SDSS optical studies of active galaxies 
( Miller^j^ilJ 200311 and bl ack hole demographics arguments 
(e.g.. iMareoni^t alJl200^) all support values of over 10 8 yr. 
In order to observe a radio galaxy when its nucleus is still 
actively feeding its jet, it must be intercepted by our light 
cone at some epoch between the cosmic time of its birth and 
the time when its beam is switched off, i.e., within an interval 
of TMaxAge after its birth. For this interception to occur the 
source must lie inside a certain cosmic volume shell, the 
"Relevant Comoving Volume Element" , Vc (BRW) . 

For a spatially flat (k = 0) universe, if r is the radial co- 
moving coordinate, Vc = 4nR 3 (t) ( rf — rfj /3, where R(t) 
is the scale factor of the universe at cosmic time t, and r\ 
and V2 are the inner and outer radial coordinates of the vol- 
ume shell (at tbirth)- The value of Vc is the relevant volume 
at the epoch Zbirth (or tbirth)- The corresponding proper vol- 
ume now is Vc(z = 0) = (1 + z birt h) 3 Vc{z birth ). 

The sources are assumed to be distributed in redshift 
according to a gaussian RLF with (BRW Eq. 24) 



p(z) oc exp 



I ( z ~ z ° 
"2 V o z 



(1) 



a distribution that peaks at zo and has standard devia- 
tion of o z ■ According to the RLF of Will ott et"afl ^ 



2.2j_a z = 0.6, and we use these values in our simula- 



tions. I^Wme^e^fl] i2004T) have given a more recent compu- 
tation of the RLF where the values are zo = 1.7, a z = 0.45 
(see their Table 5). The number of sources born at some 
cosmic time (t), per unit cosmic time, per unit comoving 
volume element at redshift zero is found from the relation 
p(t)dt — p(z)dz. For a homogeneous and isotropic universe, 
this distribution is valid at all epochs throughout the space. 
At a particular redshift Zbirth: the comoving volume (Vc) is 
found. Then multiplying Vc by p(zbi r th) gives the number of 
sources born at z b i r th (per solid angle) in the chosen inter- 
val in cosmic time which are intercepted by our light cone: 
Nbarn <x V (z = 0) p(z blrth ) . The total number, N bor n, is 
obtained by using a normalization factor in the above pro- 
portionality which takes into account the sky area of the 
observed data sample. 

Homogeneity of the universe implies that the sources 
are randomly distributed within the comoving volume shell. 
The age of a source, T age , is the time after tbirth it is inter- 
cepted by our light cone; in our computations it is drawn 
randomly from to TMaxAge, but weighted so that sources 
are distributed uniformly in volume within the comoving 
volume shell. 

In each simulation (run) we have generated a very large 
number of sources, over a wide range of cosmic time. We 
find the number of sources born at some Zbirth which will 
intercept our light cone, the age T age (denoted by t hence- 
forth) of each source, and the redshift at which we observe 
it (denoted by z henceforth), which is derived from T b 3 , the 
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cosmic time at which the light we see was emitted from the 
source. 

As very powerful sources are much rarer than weaker 
ones, each of the sources generated is assigned a jet power 
Qo (which is assumed to remain constant throughout its age) 
according to the probability distribution (BRW Eq. 38) 

p(Q())dQo OC Qq X dQ if Qmin < Qo < Qmax, 

= if Qo > Qmax or Q < Q m i„. (2) 

Here the power index x is positive, and we initially adopted 
the values used by BRW: x = 2.6, Qmin = 5 x 10 3r W, and 
Qmax — 5 x 10 42 W. Our best fit values of x are higher and 
are discussed in §5. 

An initial Monte Carlo population generation is com- 
pleted when t, z and Qo are randomly assigned to each 
source of the population according to the above prescrip- 
tions. Each source in that population is then allowed to 
evolve according to a model described in the following sec- 
tion, giving the observable quantities other than z: P, D and 
a. 



Table 1. Default Values of the Model Parameters 3, 



Parameter KDA BRW MK 



(3 1.9 1.5 1.5 



ao (kpc) 


2 


10 


10 


PO (kg m -3 ) 


7.2 x icr 22 


1.67 X 10~ 23 


1.7 x 10 


r* 


5/3 


5/3 


5/3 


r c 


4/3 


4/3 




Pb 


4/3 








1.3 






Irnin(hs) 


1 


1 


10 


^fm ax(hs) 


Infinity 


10 14 


10 7 


p 


2.14 


2.14 


2.23 


r hs (kpc) 




2.5 


2.5 


tbs (yr) 




10 s 




tbf (yr) 




1 




V 






0.1 


e 






1.0 


T 






2 x 10" 



3 MODELS OF RADIO LOBE EVOLUTION 

A sta ndard basic m o del of FR II extragalacti c radio sources 
(e.g., IScheuerlll97l iBlandford fc Reeslll974l) is widely ac- 
cepted. A powerful RG consists of the central active nu- 
cleus and two jets emerging from opposite sides of it. After 
traveling substantial distances the plasma in these jets col- 
lides with a tenuous environment. There the jets terminate 
in a shock where relativistic electrons are accelerated and 
hotspots are formed; the plasma passing through the ter- 
minal shocks inflate the huge lobes of energetic particles. A 
bow shock propagates into the surrounding gas ahead of the 
jets. 

The radio power evolution models that we compare are 
those given by KDA, BRW and MK. In brief, the physics of 
these models differ mainly in the ways in which particles are 
assumed to be transported from the jet through the hotspot 
and into the lobe. KDA assume a constant injection index, p, 
for the energy number distribution so N(E) oc E~ p , for the 
radiating relativistic particles while the particles are injected 
from the hotspots into the lobes. BRW assume that the in- 
jection index varies between the different energy regimes, 
as governed by the break frequencies discussed below. MK 
assume a constant injection index but also argue that the 
particles are re-accelerated by some turbulent process in the 
head during transport to the lobes. Several key points of 
each model and additional differences are noted in §§3.2 - 
3.4, although the reader should refer to the original papers 
for a thorough understanding of each model's details. Ta- 
ble 1 lists the default values of the major model parameters 
(those used by the authors). We varied these parameters in 
our extensive simulations described in §4.3. The only pa- 
rameter whose variation was not considered is the adiabatic 
index of the external environment, which was adopted as 
usual as T x = 5/3. 

3.1 Dynamical Expansion and Emission 

In all of the models we consider here the ambient gas around 
the double radio sources, into which the lobes propagate, is 



a See text and the original papers for parameter definitions. 

taken to have a power-law radial density distribution scaling 
with dist ance r > ap from the ce nter of the host galaxy 
(Eq. 2 of lKaiser fc Alexanderlll997tl . 



p(r) = pa 



(- 

\ao 



(3) 



where the central density, po, scale length, ao, and radial 
density index, f3, are given by the particular model. We 
follow BRW and assume that the external density profile 
is invariant with redshift. While such a typical radial den- 
sity distribution is appropriate on average for small z, this 
may not to be a good approximation at the redshifts cor- 
responding to the quasar era, which witnessed a 10 2 -10 3 
times higher co-moving density of powerful radio-loud el- 
lipticals fe.g.. ljackson fc Wallll99^ . We note that for very 
large sources the density will depart from a single power-law 
with radius and eventually approach a constant value ap- 
propriate to the intergalactic medium at that redshift (e .g., 
iGopal-Krishna fc Wiitalll987t iFurlanetto fc Loeb1l200lD . If 
the ambient density approaches a constant value at radial 
length scales of ~ 100 kpc, then radio sources grow to some- 
what smaller sizes and have larger lobe powers. We will con- 
sider this more complicated situation i n future work. 

F rom d imensional arguments llKaiser fc Alexander! 
Il997l or KA: lKomissarov fc Falie]|l998l) the total linear size 
(from one hotspot to the other) of a radio source at an age 
t can be expressed as 



D(t) = 2cia | — 

og po 



t 3 Q 



1/(5-/3) 



(4) 



where ci ~ 1, is model dependent, but weakly varying, as 
discussed below. The jump conditions at the external bow 
shock and the expression for linear size give the pressure of 
the head plasma immediately downstream of the bow shock 
as (KA Eq. 12) 



Ph(t) 



18c 2 " 5 



(r„ + 1) (5 - py 



3 3/3^,2- 
Po a o V|) 



0-, 1/(5-/3) 



(5) 
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An ensemble of 71(7) relativistic electrons with Lorentz 
factor 7 in a volume V with magnetic field B emits syn- 
chrotron power per unit frequency, per unit solid angle given 
by (KDA Eq. 2) 

P^^f'-n^V (6) 

in units of W Hz -1 sr _1 , with or the Thomson cross-section 
and po the permeability of free space. These relativistic elec- 
trons are injected into the lobe from the hotspot via the 
head, an extended region of turbulent acceleration around 
the hotspot. 



3.2 The KDA Model 

For the density profile of the external atmosphere this model 
uses po = 7.2 x 10" 22 kg m ~ 3 , a = 2 kpc and (3 - 1.9. 
These values are argued to be typi cal for an elliptical galaxy 
out to » 100 kpc from its cente r ifFonnamJon es fc Tucked 
ll985l : ICanizares. Fabbiano fc Trinchierill987l) . The factor ci 
(Eq. 4) is given by Eq. (32) of KA, which by their Eqs. (37) 
and (38) depends weakly on Rt, the axial ratio, defined as 
the ratio of the length of the source and its width. For Rt = 
1.3, the value adopted by the authors and us, ci = 1.23. 

The ratio of pressure in the head to that in the cocoon 
was taken by KDA to be (Eq. 38 of KA) p h /p c = 4R T . We 
follow thi s prescription; howev e r, the hydrodynamical simu- 
lations of iKaiser fc Alexander! <ll999T) found th is ratio to be 
an overestimate. The "improved" KDA model jKaiserll200d) 
obtai ns an empirical formula for this ratio as (Eg. 7of lKaiser1 
2000) p h / Pc = (2.14 - 0.52/3) fl 2 04-0 ' 25 ' 3 . We are exploring 
this alternative approach and will give results using it in the 
next paper in this series. 

The electrons are assumed to be accelerated in the 
hotspot at time ti, with corresponding initial Lorentz factor 
7i. The energy distribution of the electrons injected into the 
lobe is a power law function of 74, n(7i) oc % v \ p is taken to 
be constant. The electron energies evolve in time according 
to (Eq. 4 of KDA) 



(h, 



O17 4<T T -i , 

6t 3m e c 



(7) 



Here the lobe electrons undergo energy losses via adiabatic 
expansion (V oc t ai with ai = (4 + P) / [T c (5 - /3)] and 
r c the adiabatic index in the cocoon, KDA), IC scatter- 
ing off the CMB photons and synchrotron losses. The mag- 
netic field (assumed to be completely tangled) with en- 
ergy density ub and adiabatic index Fb = 4/3, satisfies 
u B oc B 2 (t) oc t~ TBai . The energy density of the CMB, u e , 
is taken to be constant for an individual radio source as each 
source evolves for only a few times 10 8 years. 

The KDA model does not distinguish between the head 
and hotspot, and considers self-similar expansion of the 
head, where the jet terminates. The cocoon is split into 
many small volume elements, each of which is allowed to 
evolve by expanding adiabatically (changing the pressure 
from head pressure ph{pi) to cocoon pressure p c {U)) and 
undergoing the various loss processes. The energy of each 
volume element in the lobe is equated to the energy it had 
while in the head minus the work done by the volume in 
adiabatically expanding from the head to the lobe. The ra- 
dio emission from such a volume element is calculated, using 



the expressions of cocoon pressure and the energy distribu- 
tion function. The total emission at a frequency v is then 
obtained by summing over the contributions from all such 
small elements in the lobe. The expression for P v , given by 
Eq. (16) of KDA is a complicated integration over injection 
time ti. This integration being analytically intractable, we 
used numerical techniques to get P v . 



3.3 The BRW model 



The ambient gas 



are po 
P = 1.5. 



density 
-23 



parameters adopted by BRW 
1.67 x 10"* 3 kg m -3 , a = 10 kpc and 
These are based o n polarization measurement s 



of lobe synchrotron emission iGarrington fc C onwaj|ll9^|) 
and X-ray images of massiv e ellipticals (e.g.. ISarazi nl ll988l : 
iMulchaev fc Zabludofdll998ft . A value of a = 1.8 is adopted 
in Eqs. (4) and (5), as BRW found it to give the best fit 
between models and data. 

This model assumes the hotspot to be a compact region 
(the working surface moving around as in Scheuer's ( 1982) 
"dentist's drill model") within the whole head region. Con- 
sidering the expansion o f the head and its bow shock (also 
iBegelman fc Cioffilll989l) . the environmental ram pressure is 
related to the average internal pressure in the head (Eq. 5). 
The pressure in the lobe in taken to be a constant factor 
(1/6) of the head pressure. 

The jet, of constant bulk power Qo, terminates at the 
hotspot, which is taken to be of constant radius, r^ s = 2.5 
kpc, in BRW. The pressure in the hotspot, phs, is given 
by the stagnation pressure in the post-jet shock, pn s = 
Qo/icAhs)- Here A^ s {= ""'"/L) i s the area normal to the 
jet over which the jet thrust operates. The hotspot mag- 
netic field, assumed to be tangled, is given by B^ s = 
3/J.oQo/(cAh s ), where the equipartition assumption has been 
made. The break frequency for synchrotron radiation in the 
hotspot is (Eq. 12 of BRW) 



Vbh = 



9c 7 B ha 



4(B£, 



4- R 2 

T D CMB 



(8) 



where c 7 is 1.12 x 10 3 nT 3 Myr 2 GHz l)Leahvlll99lL and 
the equivalent magnetic field due to the CMB is Bcmb = 
0.318(1 + z) 2 nT. The synchrotron age, t s , of the electron 
population is determined by the length of their exposure to 
the hotspot magnetic field before they reach the lobe. This 
longest dwell time in the hotspot is taken as tb a = 10 s yr, 
and the shortest dwell time tbf = 1 yr. 

In §8.4.2 of BRW it is shown that this model roughly 
follows the KDA prescription of lobe luminosity, but with 
two main differences. First, while the particles are injected 
from the hotspot to the lobe, the injection index is governed 
by the breaks in the energy distribution of particles (unlike 
the constant injection index of KDA). Second, the constant 
hotspot pressure governs the adiabatic expansion losses out 
of the hotspot (for particles injected into the lobe), while 
in KDA the head pressure (which evolves with time) drove 
the adiabatic losses. In BRW the head pressure (Eq. 5) only 
drives the source expansion. 

The details of the energy distribution (as a function of 
the Lorentz factor) of particles in the hotspot are shown in 
Fig. 11 of BRW. Our calculations usually are done assuming 
the minimum and maximum values of the particle Lorentz 



factors in the hotspot quoted by BRW: 7 



min(hs) 



1 and 
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Imax(hs) = 10 , although we have examined variations in 
these parameters. 

A population in the lobe which emits at a time T b B 
(when it intercepts our light cone), consists of particles in- 
jected from the hotspot between a time t m in (those with 
largest Lorentz factors) and T t s (smallest Lorentz factors). 
The time t m in (found following the prescription in KDA) is 
the minimum time of injection (found for every T b s ), when 
particles can still contribute to the radiation at v. 

The final expression for the power emitted (P v ) by a 
radio source at a frequency v is given by the complicated 
Eq. (21) of BRW, which we will not reproduce here. We 
solved this equation numerically. 

3.4 The MK model 

The lManolakou fc Kirkl (2002) paper employs the same ex- 
ternal density profile and source linear size expansion as does 
BRW. 

The MK model essentially follows the common prescrip- 
tions of KDA and BRW for lobe luminosity evolution, with 
the key difference involving the particle transport mecha- 
nism. Two cases are considered for the propagation of par- 
ticles from the termination shock through the hotspot and 
into the cocoon. In MK's Case A, the whole adiabatic loss 
between the hotspot and lobe (due to the pressure differ- 
ence) is computed. However, the authors found that this 
produced P — D tracks which conflicted with the observa- 
tional data. So they considered a Case B, which involves 
some re-acceleration process in the turbulent head region, 
whereby the adiabatic losses are partially compensated; MK 
found such a model is a qualitatively better fit to the data. 
Thus we consider only the case B (with re-acceleration) of 
the MK model in our present work. 

This model assumes that electrons are accelerated by 
the first-order Fermi mechanism at the jet termination shock 
and are injected into the plasma behind the shock following 
a power-law energy distribution with a constant injection 
index p. A fraction r\ of the jet power is assumed to be 
transferred into the accelerated particles at the termination 
shock. If the bul k Lorentz factor of the jet, 7j e i ~ 10, then 
2 < p < 2.3 ('e.g.. lAchterberg et al .1200 J) . The correct upper 
and lower limits of particle Lorentz factors "f m in and 'jmax 
are not obvious; MK adopt 7 m i n = Jjet- The authors say 
the results are not sensitive to j ma x', however, our different 
conclusions on this point are given in §§5-6. 

After being dumped in the primary hotspot by the jet, 
the electrons encounter turbulent motions of the plasma in 
transit through the head and finally reach the lobe. In this 
transition through the head the electrons are subject to syn- 
chrotron losses (in the strong magnetic field behind the ter- 
mination shock) and IC losses off the CMB. The effects of 
the losses depend on the distribution of the "escape times" , 
i.e., the probability distribution of how many particles es- 
cape after a certain time interval. A generalized transport 
process is considered, with e (denoted as a in MK) being 
the transport parameter (or the diffusion index). The mean 
square distance traveled by a particle, (Ar 2 ) tx t e , with 
< e < 2. In the standard diffusive case, e = 1, with sub- 
(supra-) diffusive cases being e < 1 (> 1). Another new pa- 
rameter in this model is r, the ratio of the diffusive transport 
time and cooling time of a particle at 7 m i„. 



During the transport of particles from hotspots to lobes, 
the details of re-acceleration by various proces ses have been 
considered by many previous authors (e.g. | Spruit| [l98q: 
Bceclman & Kirk 1990; Mano lakou. Anastasiadis fc Vlahosl 
Il999l iGieseler fc Jonesll2000l) . MK simply assume that in 
the presence of reacceleration, the distribution of electrons 
entering the lobe is described by a power law above a lower 
cut-off energy, and at higher energies is modified by syn- 
chrotron and IC losses. 

Once the electrons have reached the lobe, they undergo 
adiabatic, IC and synchrotron losses, similarly to the other 
models, and their energy evolution is given by Eq. 22 of 
MK. Similarly to the KDA model, for every time instant 
t, when radiating particles have Lorentz factor 7, there is 
an earliest time, ti, at which particles injected into the lobe 
contribute to the radiation at t; ti can be obtained following 
the prescription given in Eqs. (25) and (26) of MK. 

The final expression for power emitted at a frequency v, 
P v is given in Eq. (27) of MK. The authors used r^s ~ 2.5 
kpc as the hotspot radius; however, we find that the MK 
model power results are actually independent of the hotspot 
area Ah s . From Eq. (2) in MK, Uh oc and from their 

Eq. (6), t oc A^/f", where a = (4 + /?) / (5 — Hence ui obe 
(MK Eq. 5), b s (MK Eq. 22) and P v are independent of A hs . 



4 OBSERVATIONAL AND SIMULATED 
SAMPLES 

4.1 Selection Criteria and Observed 
Characteristics 

These models predict the emission from the radio lobes, 
which are taken to (and usually do) dominate the emission 
from extended FR II RGs. As is well known and is discussed 
in detail in BRW, at relatively low frequencies (~ 151 MHz) 
the radio flux observed is predominantly the emission from 
the cocoon or the lobe (with negligible contribution from the 
hotspots, jets or nucleus), and so these evolutionary mod- 
els should fit the data best at such frequencies. At GHz 
frequencies, substantial contributions from Doppler boosted 
core or jet emission would often be present, especially for 
old quasars, but the slowly advancing lobes will still emit 
nearly isotropically. In addition, at these higher frequencies 
the effects of synchrotron, adiabatic and IC losses are more 
severe. At very low frequencies (< 100 MHz), there are extra 
complications affecting the emission from synchrotron self- 
absorption, free-free absorption, and the poorly known low 
energy cut-off to the relativistic synchrotron emitting parti- 
cles. Therefore samples such as those produced by the Cam- 
bridge group over the past decades, which were observed at 
between 151 and 178 MHz, and cover much of the northern 
sky, are most appropriate for this work. 

We adopt observational samples from complete radio 
surveys (Table 2), each of which contains all the radio 
sources within each survey's flux limits and which are found 
inside smaller sky areas, for deeper surveys. Redshifts have 
been obtained for the great majority of these radio sources. 
This lower flux limit brings in a P — z correlation, since P de- 
creases as z increases. To decouple this P — z correlation one 
must use multiple complete samples at increasingly fainter 
flux limits. 
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Table 2. Observational Samples 



Survey 


Flux Limit 

(Jy) 


No. of Sources a 


Sky Area 
(sr) 


3CRR 


5i78 b > 10.9 


145 


4.23 




Sisi > 12.4 






6CE 


2 5isi ^ 3.93 


56 


0.102 


7CRS 


Sisi > 0.5 


126 


0.022 


7CI 


Sxsx > 0.51 


37 


0.0061 


7CII 


Sisi ^ 0.48 


37 


0.0069 


7CIII 


5X51 > 0.5 


52 


0.009 



a Only FR II RGs considered. 

'Flux at 178 MHz, the frequency at which the 3CRR survey was 
performed. Sns for these sources are converted to flux at 151 
MHz, S151, using a constant spectral index of 0.8. 

For an individual source in each survey, the following 
characteristics were considered: the redshift (z), the specific 
power at 151 MHz (Pxsx) in W Hz -1 sr -1 , the total pro- 
jected linear size (D) in kpc, and the spectral index at 151 
MHz (0151) converted to the rest frame of the source. The 
redshifts in the samples are spectroscopically determined for 
the vast majority of the sources. For the 3CRR catalog the 
redshift completeness is 100%, for 6CE it is 98% and it is 
92% for 7CRS. 



4.2 Observational Sample Details 

Henceforth, 3C, 6C, and 7C refer to the refined surveys 
3CRR, 6CE and 7CRS, respectively, as described below. We 
excluded FR I RGs from the following catalogs and consid- 
ered only FR II sources (including quasars, weak quasars, 
high-excitation FR II RGs and low-excitation FR II RGs). 

3CRR: This is the Third Cambridge Re- 
vised Revised sample of extragalactic radio sources 
jLaing. Rilev fc Longairlll983lh We adopted the data from 
the online compilation of the list by Willott 1 . In 3CRR 
the observations were done at a frequency of 178 MHz, so 
for each 3CRR source P178 (specific power at 178 MHz) 
was obtained and then converted to Pxsx using a standard 
average spectral index of 0.8. Given the closeness of these 
two frequencies, any reasonable variations in a would make 
for only small differences in the derived P151 v alues. 

6CE: The Sixth Cambridge radio survey bv lEalesI Jl985T) 
is the original 6C survey. We ado pt the sample from the 
reselec ted and updated version in iRawlings. Eales fc Lacvl 
teOOlf) . along with the most recent redshifts, which have 
been updated online by Rawlings 2 . 

7CRS: The Seventh Cambridge Redshift Survey is a 
combination of parts I, II and III of the original 7C survey 
jMcGilchrist et alll99o]). F or 7C-I and II we adopt Pi 5 i and 
z from Wi llott et alJ a2003l) : (their Tables 2 and 3, which use 



1 http:/ /www-astro. physics. ox. ac.uk/~cjw/3crr/3crr. html 

2 http:/ /www-astro. physics. ox. ac.uk/~sr/6ce. html 



the present consensus cosmology). The values of D were ob- 
tained from a web-site maintained by Steve Rawlings 3 ; how- 
ever, the axsi values are not available in a collated form in 
the literature, and only a few individual sources have these 
values published. Thus we used CV151 for 7C-III only. For 
7C-III, the reduced data, including redshift, flux density in 
Jy, angular size in arcsec and spectral index between 38 and 
151 MHz were kindly provided to us by Chris Willott; from 
these we computed the relevant observational parameters in 
the cosmology we use. The observed spectral index between 
38 and 151 MHz was taken as the rest-frame 151 MHz spec- 
tral index (a fairly good estimate, at least for the higher 
z sources). The relevant sample can be found in Table 9 
(conta ining both the 7CIII and NEC samples) of lLacv et alJ 
( 1999) or online from the website of Oxford University 4 (but 
with a different cosmology). 

4.3 The Simulated Surveys 

Large radio source populations are randomly generated, ac- 
cording to the source age, redshift and beam power distri- 
butions as given in §2, for each choice of model parameters. 
Each simulated source, in a population, is then allowed to 
evolve in age according to a power evolution model discussed 
in §3. The evolution is to be done at a rest frame frequency 
of the source, so if the frequency of observation is v t s = 151 
MHz, and a source is observed at redshift z, then it is evolved 
at a frequency v res t = 151 x (1 + z) MHz. 

The monochromatic power (Pxsx in W Hz -1 sr -1 ) each 
source would emit at the T (, s corresponding to it is calcu- 
lated; this depends on the model as described in §3. At this 
cosmic time (T ba)% its redshift, and hence its distance from 
us, is found. The flux (in units of Jy = 10~ 26 W Hz" 1 m~ 2 ) 
of this source is then obtained fusing IPeacockl jl999f) : Eqs. 
3.87, 3.76 and 3.10 for a flat universe), given that it emitted 
Pxsi from the cosmic distance calculated. If the flux for a 
source is greater than a (lower) survey flux limit (or between 
two flux limits in the case of 6C) then that source is consid- 
ered to be detected in the corresponding simulated survey, 
and counted for the later comparisons with real data. 

It is assumed in our simulations that the radio jets feed- 
ing the lobes (or the central AGN) stay "on" only for the 
time T ag e corresponding to each source (§2), which is also 
taken as the lifetime of a source. After the time T age , the 
relativistic plasma in the lobes continue to radiate but the 
flux drops very rapidly once the central engine has stopped 
to feed the lobes. So the sources can be considered to be 
turned "off" instantaneously after T age . This assumption is 
supported by the fact that the radio powers (P151) drop sub- 
stantially while the jets are still on (i.e., within the time T age 
after birth), as shown by the P — D tracks in §5.1. 

To perform our simulations we initially generate an en- 
semble of a huge number (a few 10 6 ) of pseudo-radio sources. 
After evolving each source by the above prescription, the en- 
semble is examined to see how many of them would actually 
be detected in a simulated complete survey. The population 
size is then chosen for this parameter set in the next run so 
as to get a comparable number detected in the simulations 

3 http:/ /www-astro. physics. ox. ac.uk/~sr/grimestable.ascii 

4 http:/ /www-astro. physics. ox. ac.uk/~cjw/7crs/7crs. html 
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as are found in the real surveys. To do this, we usually had 
to generate such "standard" ensembles containing ~ 10 6 to 
10 7 radio sources. Assuming the observed regions are fair 
samples of the universe, the population size is proportional 
to the sky area of a survey. The populations needed in or- 
der to simulate the 6C and 7C surveys are generated from 
that of 3C by reducing the total 3C population size accord- 
ing to the corresponding sky area ratio. Given a 3C initial 
ensemble of size S3C, the populations for 6C and 7C are 
created by plucking sources randomly from that initial en- 
semble, and producing populations of sizes Sec = S , 3c/41.5 
and SVc = 5 , 3 c/192.3. 

The initial populations generated for comparison with 
the 6C and 7C data following the above procedure detected 
more or less comparable numbers of sources in the simula- 
tions as compared to the actual surveys. We compute the 
over (-under) detection factors, defined as the ratios of the 
number of sources detected in the simulated 6C and 7C sur- 
veys to the numbers in the actual surveys divided by the 
same ratio for the 3C survey. The deviation of these ratios 
from 1.0 (see discussion in §6) may be considered a measure- 
ment of the statistical (sample) variance. 

Each model gives the total linear size of a radio source, 
but we observe each one as projected on the plane of the sky. 
This is incorporated into the simulations as follows. Sources 
are considered to be randomly oriented in the sky, with the 
angle to the line of sight (8) of each source drawn from 
a distribution uniform in (1 — cos 8). The projected length 
of each simulated source is then, D pro j = D(t) x sinf? = 
D(t) x -y/rjv(2 — rjv), where rjv is an uniform random num- 
ber between and 1, and D(t) is the total linear size of the 
source. For compactness, hereafter we denote the projected 
size Dp ro j as just D. 

The initial population of sources was generated and 
each lobe power evolution model was implemented in C, 
and the other suppor ting codes were w ritten in IDL. Nu- 
merical Recipes in C jPress et alj l2002) were used to speed 
up the calculations of lobe powers for the huge ensembles of 
sources. 

In doing the statistical tests (§5.2 and §5.3) we com- 
pared the model predictions with the observational samples 
as follows. In a single run a random initial population of 
millions of sources was generated such that after evolution 
of each source in the ensemble and after comparing each 
to the flux limits, the ensemble produced simulated sam- 
ples (for the 3C, 6C and 7C catalogs) which were of sizes 
comparable to or larger than the real surveys. The simulated 
samples were then reduced in size, if necessary, by uniformly 
selecting sources from them. In particular, this was done by 
selecting every (A r s i m /A r samp )'th source from a simulated 
survey, where N 3arnp is the number of sources in one of the 
real surveys 3C, 6C or 7C and N sim (usually > N sa m P ), is 
the number of sources in the simulated survey. Finally sta- 
tistical tests (whose results are tabulated) were done on the 
[P, D, z, a] data from the real surveys and a similar sized 
simulated sample generated from a single random seed. 

Each of the [P — D — z — a] plane figures, described 
in §6, show the final simulated sample (after reduction to 
the actual data sample sizes) of the random run done us- 
ing specific parameters for each semi-analytical model. The 
plotted cases were among the best overall fits for each model 
as determined by the KS tests. 
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Figure 1. P — D tracks of three sources with jet powers (Qo) 
in Watts and redshifts (z) of [1.3 X 10 40 ,2.0], [1.3 X 10 39 ,0.5], 
[1.3 X 10 38 , 0.2] (from top to bottom). Each of the dashed, dotted 
and solid curves correspond to the tracks predicted by the default 
versions of the BRW, KDA and MK models respectively. The 
crosses on the tracks denote source lifetimes of 1, 10, 20, 30, 
90, 100, 150, 200, 250, 300 Myr. 

5 RESULTS 

5.1 P-D Tracks 

As a radio source gets older, its power (P) vs. linear-size 
(D) track becomes steeper. While this is true for all mod- 
els, the rate of steepening is different in the three models, 
as seen from Fig. 1. These P-D tracks have been generated 
using the default parameters of each model (given in Table 
1), by allowing each source (with beam powers and redshifts 
given in the plot) to evolve at frequency v — 151 MHz. 
For this Figure (alone) the total linear sizes were converted 
to the projected sizes assuming an average viewing angle 
to the line of sight of 39.5° (following KDA). These tracks 
are in agreement with the conclusion drawn by MK that 
their P-D tracks are more akin to those presented by KDA. 
Crude evaluations of the quality of different models and the 
allowable ranges of parameters for them can be found by 
comparing the regions in the P-D diagram that are actu- 
ally populated with those that are accessible to models with 
those parameters (e.g., KDA, MK). 

On examining different tracks, it is found that the lu- 
minosity falls off faster for sources with high power and 
high redshift. The higher the redshift of a radio source, 
the shorter the fraction of its life which can be detected 
by flux limited radio surveys. This p oint was noted by 
Gopal-Krishna._Wiita_fc Sari rjallil dl989T) and was stressed 
bv lBluru1^ii^^R^wlrngsri^99l^ who coined the phrase 
youth-redshift degeneracy to describe it. 

5.2 1-Dimensional Kolmogorov-Smirnov Tests 

5.2.1 Statistical Tests and Default Parameter Results 

For our first attempt to quantitatively compare the simu- 
lated radio surveys to the actual data, we used 1-dimensional 
Kolmogorov-Smirnov (1-D KS) statistical tests. Based on 
the results of such tests we chose some parameter variations 
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for the models on which two more statistical tests were done 
(§5.3). 

In the 1-D KS tests, each of the distribution's key char- 
acteristics [P, D, z, a] of the radio sources detected in the 
simulated surveys were compared to those of the sources in 
the real radio surveys 3C, 6C, and 7C, according to proce- 
dures given at the end of §4.3. The KS probabilities, V, that 
the two data sets being compared are drawn from the same 
distribution function, were taken to be a figure of merit of 
each model used in the simulation. High values of V (close 
to 1.0) indicate good fit, and very small V imply that the 
model and data distributions are significantly different. We 
consider twelve test statistics in total (the twelve probabili- 
ties found from the KS statistics for comparisons of each of 
P, D, z, a for each of the three radio surveys), which quan- 
tifies the closeness of the model fits to the data. 

In order to quantify the overall success of a model we 
would prefer to have a single figure-of-merit instead of twelve 
individual ones, but there is no obvious way to produce such 
a statistic, particularly since the three surveys have signif- 
icantly different numbers of objects. A likelihood type test 
would involve the product instead of the sum of the KS 
probabilities, but given the extremely small values of these 
products we rejected this figure of merit as not providing 
a useful discrimination between the models. Here we have 
combined the 1-D KS probabilities in two ways. 

First, we add the KS probability statistic for compar- 
isons of P,D,z,a (i.e. P(P) + P(D) + V(z) + V[a)) for 
the three surveys, weighting the statistic of a survey by the 
square-root of the number of simulated sources detected in 
that survey. So the first overall figure of merit of a model, 
which we denote as r P[p 1 D,z,a], is given by: 

TlP,D,„ a] = [P(P) + V(D) + V(z) + T(a)] 3C + 
[V(P) + V{D) + V{z) + V(a)] 6C + 

^^[V(P)+V(D)+r(z)+r(a)] 7C , (9) 

where N3C, Nqc and Nye are, respectively, the number of 
sources detected in each of the simulated surveys with a 
particular parameter set used in the model. As noted above, 
if the simulations "detect" too many sources as compared 
to the data, then each of the resulting simulation survey 
samples for 3C, 6C, 7C are reduced by uniformly removing 
sources to make the final simulation sample sizes equal to 
that of the data samples. 

The second figure of merit we employ adds the KS 
statistic probabilities for P and z to twice the probability 
for D, i.e. V(P) + 2V(D) + P(z) + V(a) for the three sur- 
veys, using the same weighting method. We denote this as 
V[p^2d,z,c\- This second choice was considered because re- 
sults for P and z usually correlate (due to flux-limit argu- 
ments); thus double weighting the probability for D dilutes 
the impact of the [P(P),P(z)} correlation. Unsurprisingly, 
in most of the runs we have performed the combined test 
statistics V[p,D,x,a] and V\p^D,z,a\ behaved in a similar fash- 
ion. 

Unfortunately, for the complicated problem of radio 
source cosmological evolution, which involves many param- 
eters and several dimensions, any figure of merit based upon 
1-D KS tests is a crude approach in comparing models 



with observations. We attempted to use multi-dimensional 
statistical tests (e.g.. [Holmstrom. Sain fc Miettinenl ll99Et 
lLoudin fc Miettinenl I2003T) which, in principle, could yield 
a more robust single figure of merit for the fit of our dis- 
tributions to the data in [P, D, z, a] space. Unfortunately, 
the limited sizes of the observational samples (< 150) pre- 
clude obtaining reliable results from such generalizations of 
the KS test. Here we are trying to fit four variables, namely 
[P, D, z,a]; in practice the minimum useful sample size re- 
quired would be ~ 10 4 for a four- dimensional test. In future 
work we plan to expand our method to include simulations 
of large scale radio surveys containing many thousands of 
sources, such as NVSS and FIRST, which can be made ad- 
equate ly complete in z t hrough optical identifications from 
SDSS ( Ivezic et alJl2004) . Then we might successfully incor- 
porate a multi-dimensional test. 

We call the parameters governing the power evolution 
as given in the papers KDA, BRW (including the parameters 
used for the initial radio source population generation) and 
MK to be the default ones. For KDA and MK the authors 
discuss some alternate parameter sets in the respective pa- 
pers, our defaults are their first and favored parameters. To 
save space, the statistics of only a subset of all the runs we 
have performed are shown in the present work. We include 
the cases which give the best 1-D KS statistical results. 

We present the KS test results in tables grouped by 
radio source evolution model, with each entry illustrating 
a different parameter set. Table 3 gives our results for the 
KDA model, Table 4 for the BRW model, and Table 5 for 
the MK model. The tables for each model follow the same 
format and pattern. Hence we describe only the table entries 
for the KDA model (Table 3). 

Each of the Tables 3, 4 and 5 give the individual KS 
statistic probabilities V(P), V(D), V(z) and V(a) for some 
of the initial runs. The results for each model are given in 
three consecutive rows. The first column lists the values of 
the RG source distribution index x and TuaxAge (in Myr) 
used for the initial population generation in that model run; 
these two parameters were expected to be most important 
in governing the numbers of acceptable sources each Monte 
Carlo simulation would generate. The second column first 
lists the parameter(s) which has (have) been varied from 
the default case in the top row(s), and then gives the initial 
population (ensemble) size used for 3C simulation in that 
model. The third column notes to which survey KS proba- 
bilities given in the next columns correspond. The first row 
in columns 3, 4, 5, 6 and 7 gives the 3C results, with the 
second and third rows giving the 6C and 7C results, respec- 
tively. The fourth, fifth, sixth and seventh columns show the 
values of V(P), V(D), V(z) and V(a) respectively for each 
of the surveys, 3C, 6C and 7C. The final, eighth column, 
lists the combined probabilities, "Pyp^D,?,^ and P[j>2.D,z,a]> 
in two consecutive rows, for each particular parameter set. 

To begin with, an initial population, generated using 
the default parameters from BRW for RG population gen- 
eration, was evolved according to the three different default 
models discussed before. The simulated sources detected (ac- 
cording to the prescription in §4.3) were compared to the ac- 
tual data in the 3C, 6C, and 7C catalogs. As shown by the 
KS test statistics of the first three rows of tables 3, 4 and 5, 
the model fits are all very poor. The main problem is that 
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too many high- 2 and too few low-z sources were produced 
by the models as compared to the data. 

5.2.2 Dependence on Source Slope Parameter, x 

To look for improved agreement between simulation and 
data, we decided to steepen the beam power distribution 
function of the sources generated in the initial population. 
This significant modification was expected to produce fewer 
high P - high z sources, and the exponent in the power law 
distribution of the jet powers, x, was increased from x — 2.6 
(as used by BRW) in intervals of Ax = 0.2 or 0.3. For the 
KDA and MK models the overall statistics improved the 
most at x = 3.0, but were less good for x = 3.3 or 3.6. For 
BRW the P and z fits were best for x = 3.6, making the 
overall performance look fairly good, but the D fits were 
all very bad. As will be discussed further below, the former 
modification (a; = 3.0) produced a clear overall improvement 
for the BRW model too. 

The initial population generated with x = 3 (but oth- 
erwise using the BRW prescription), was evolved accord- 
ing to the KDA and MK models. The corresponding KS 
statistics are given as the third entries in Tables 3, 4 and 
5. For the BRW model, the big population generated using 
x — 3.6, had a very strong P — D anticorrelation, producing 
too many small sources and too few large ones. The com- 
bined KS statistics were also much worse than those of the 
KDA and MK models, so we do not list any BRW model 
results with x > 3.0. 

Some of the 12 KS probabilities for the KDA and MK 
models (albeit very few for BRW) provide acceptable fits to 
the data. To search for possible further improvements we 
varied the other parameters describing the power evolution 
in the models as described below. 

Accepting x = 3.0 as a tentative value for the exponent 
of the beam power distribution for the generated initial pop- 
ulation, we then varied the parameters governing the lobe 
power evolution of KDA and MK models. For BRW the ex- 
ponent x — 3.6 was initially accepted as it gave good fits 
for P and z, though as noted above, the D fit was very 
poor, and we do not display these results. Simulations were 
done by setting the parameter values at the end points of 
physically reasonable ranges; for example we might perform 
two additional runs using the same initial population but we 
would set a parameter to half or twice its default value. 

Simulated surveys were constructed using the param- 
eter listing given in Table 3 (each variation done one at a 
time) for the KDA power evolution model. Simulations done 
with higher axial ratios (Rt = 2.0, 2.5, 3.0, 4.0, 5.0) which 
are favored by morphological data, all yielded severe under- 
detections when compared to the actual number of sources 
in the catalog. Hence we adopted the default value of the 
axial ratio, Rt = 1.3, as did KDA. 

The results for the changes of parameters considered for 
the MK power evolution model are given in Table 5. 

As seen from the tables, several of the 12 KS proba- 
bilities for some cases give acceptable fits, but it is diffi- 
cult to find a single model where all are really good fits. In 
other words, none of the models discussed here simultane- 
ously provide good fits to the data from all of the three radio 
surveys considered. As noted above, P and z seem to cor- 
relate together in most cases because they are related when 



we pick up radio sources by imposing a flux limit on them. 
In some cases P and/or z fits are good and those to D are 
bad; and vice versa. The fits to a are almost always poor. 
The KS statistics for model runs which gave any further im- 
provement over the "improved" default case (x = 3, Default 
Model Parameters) can be found from Tables 3-5. 

5.2.3 Dependence on RG Maximum Age 

An important parameter for the generation of the initial 
population of radio sources according to the BRW prescrip- 
tion is TjMaxAge. It defines the mean active lifetime of the 
RG central engine and how long the radio lobes (being fed 
by jets powered by AGN activity) continue to expand. Hence 
it is one of the most important parameters to constrain if 
we are to estimate the fraction of the relevant volume of the 
universe occupied by radio galaxies during the quasar epoch 
(§1). As our ultimate goal involves this relevant volume frac- 
tion, we aim to find the value of TMaxAge which gives the 
best fit to the data for each of the RG evolution models. We 
performed simulation runs with default parameters for each 
of the models, using initial populations with x = 3 (which 
gives the least bad fits); we then set TMaxAge to values in 
the range 50-600 Myr (in intervals of 50 Myr), and obtained 
the following results. 

For the KDA model, the combined KS probabilities, 
T > [p : r>,z,a] or P[p,2D,z,a] lacked a single maximum over the 
range in maximum age considered, and peaked at both 150 
Myr and 500 Myr. However the higher peak was adopted, 
and hence the adopted best maximum age is TMaxAge = 150 
Myr. In the other two models, the combined KS probabili- 
ties, P[p,D,z.a] or 'P[p,2D,z,ci\, varied smoothly over the range 
in maximum age considered. In the BRW model the single 
peak was at TMaxAge = 250 Myr, and in the MK model it 
was at TMaxAge = 150 Myr; hence these were adopted for 
the subsequent runs. 

Monte Carlo runs were done with the above best 
TMaxAge for each model and with x = 2.6 (the default 
from BRW) , to check if that was better. For BRW the best 
TMaxAge when combined with x — 3.0, produced better 
statistics (less bad D fit), and was hence adopted for later 
runs. In all cases x = 3.0 was better than x = 2.6 by a 
significant margin. The supporting KS statistics are given 
in Tables 3, 4 and 5 for KDA, BRW and MK, respectively. 
Hence we used initial populations with x — 3.0 and the 
above "optimal" TMaxAge values for each model in subse- 
quent runs. 

During the simulation runs we found that a few (for 
KDA and MK models) and some more (for BRW) very small 
sources (D < 1 kpc) were being "detected" in the three 
modeled surveys (mostly in 3C). The actual survey data 
has negligible numbers of such small sources which would 
not normally be considered FR II types. As any such small 
source (with D < 1 kpc) will not be regarded as a FR II 
radio galaxy, we decided to put a linear size cut-off in our 
simulations. For the KDA and MK models a cut-off of 1 kpc 
was adopted. For the BRW model we found that a cut-off 
of 10 kpc gave much better fits than did a 1 kpc cut-off. So 
in BRW we considered sources only with total linear sizes 
greater than 10 kpc. The KDA and MK simulations did not 
produce many sources with linear size < 10 kpc, hence it 
did not make much of a difference if we imposed a 10 kpc 
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or a 1 kpc size cut-off. In the results presented henceforth, 
these D cut-offs have been incorporated. 

5.2.4 Dependences on Other Model Parameters 

In order to further explore the parameter space of the mod- 
els in search of better fits, initial ensembles were generated 
using x = 3.0, TMaxAge = 150 Myr for the KDA and MK, 
and x = 3.0, TMaxAge = 250 Myr for the BRW simulations 
(see §5.2.3). The following prespription was then followed 
for each. The sources in one large random population were 
evolved several times, according to one of the three radio 
lobe power evolution models. During each evolution one of 
the model parameter values was varied around its default 
value (as in §5.2.2). All the parameters of each model given 
in Table 1 were varied, with only one variation per evolution 
run. 

Only those parameter variations that gave any improve- 
ment in statistics over the default parameter case of the same 
model, or were essentially as good as the default, were con- 
sidered further. For these parameter sets three more initial 
populations were generated having the same size and the 
same x and TMaxAge values for the different models, but 
with different pseudo-random seeds. These additional en- 
sembles were then evolved using the "improved" parameter 
sets for each model, and the 1-D KS statistics were found. 
At this point we had the KS test results for a set of four 
simulations of each of the "improved" parameter variations. 

The three cases involving variations of a single param- 
eter (previous paragraph) which gave the best statistics 
(highest mean T\pD,»,a] °f f fle 4 runs) were then found. Sim- 
ulations were then performed in which two of those parame- 
ter changes giving better fits were simultaneously employed. 
If these "2-change" variations continued to give better per- 
formances, all three changes were incorporated together in 
a single run, to see if yet better fits could be obtained. 

5.2.5 Spectral Index (a.) Behavior 

The spectral index (a) at the rest frame of a source at 151 
MHz was estimated for each source in the simulated surveys 
by considering log [v (MHz)] as the independent variable and 
\og(P v ) as the dependent. The specific powers at the T b s 
corresponding to the source (§4.3) were calculated at three 
frequencies, namely, 151, 151/(1 + z) and 151(1 + z) MHz. 
A quadratic polynomial was fitted to the log(P„) vs. log [v 
(MHz)] data. The fit coefficients a\ and a<2 where log P v = 
ao + oi log v + a,2 (log v) 1 , were obtained. These were used to 
find the spectral index as a = —ai — 2a-2 log(151/(1.0 + z)). 

The KS tests for the fits for the spectral index (a) for 
all surveys employing each model are uniformly bad (as in- 
dicated from the V(a) values in Tables 3-5). The poor 
qualit ative fits to the a distribution were already noted by 
iBlundell et"afl il999T> for their models. Still, it is the BRW 
model which gives the least unsatisfactory KS statistics for 
a fits. 

The KS statistics for a fits were extremely bad for the 
KDA model. Here, the spectral index distributions consist of 
a dense cluster at a ~ 0.58, with no sources having smaller 
a while some have steeper spectral indices up to a ~ 1.0. 
There is a weak a — D anti-correlation until D ~ 10 3 kpc, 



after which there is a trend of increasing a as D increases; 
but this involves only a few giant sources. 

The BRW model also produced mostly very poor a fits, 
but occasionally it gave quasi-acceptable KS statistics, with 
'P(a) ~ 0.01. Here, the spectral indices are almost uniformly 
distributed within a ~ 0.58 — 0.85, with some sources at 
smaller a. There is also a weak a — D anti-correlation in 
the BRW model which extends throughout the simulated 
results. 

The MK model produced the worst KS statistics for a. 
Here, the spectral indices came out very steep, with a > 0.9 
almost always found. The distribution includes a cluster at 
a ~ 0.9 — 1.0, and an extension to very steep spectra (a ~ 
1.5). Here there is a clear trend for a to be higher as D 
increases in the simulations for all three catalogs. 

Thus, it is clear that all of the models considered to date 
require modifications if they are to produce adequate rep- 
resentations of the observed radio spectral indices. Making 
such modifications is a key goal of our future work. 



5.3 Additional Statistical Tests 

In order to check the robustness of the quantitative tests 
we performed some additional statistical analyses. We se- 
lected the cases of parameter variations that gave the highest 
combined probability, P[p.D,z,a], of each model, according to 
the amplified 1-D KS test results (described in §5.2.4). We 
compared these nominally superior parameter sets for each 
model with the default versions (those with no parameter 
changes) by performing additional statistical tests on them. 



5.3.1 2-Dimensional Kolmogorov-Smirnov Tests 

We used the 2 - dimen sional (2-D) KS test procedure 
from Press et al.| (| 2002f) . wh ich is based on the work of 
iFasano fc Franceschini| jl987|). which is a variant of an ear- 
lier idea due to|PeacockJ <I1983T) . The relevant 2-dimensional 
2-sample KS probabilities (or the significance level indicat- 
ing that the two populations are drawn from the same dis- 
tribution), V, give a quantitative measure of the model fits. 
The comparisons of the model simulated samples to the real 
data samples are done in a way analogous to that for the 
1-D KS tests (§4.3). 

The 2-D KS probabilities for comparisons of the proper- 
ties P, D, z and a, taken two at a time, for the data and the 
models, were computed. Table 6 shows results for both the 
default versions and the parameter sets giving the highest 
total 1-D KS probability, denoted as varied in the table. The 
results are listed in a similar way as are the 1-D KS statis- 
tics in previous tables. The first column gives the model 
and parameter variation (if any). The third, fourth, fifth, 
sixth, seventh and eighth columns list the KS probabilities 
for comparisons of [P — z] , [P — D], [z — D], [P — a], [z — a] 
and [D — a] respectively; in each case the three rows give 
results for 3C, 6C and 7C, respectively. 

It is non-trivial to compare the models as there are 18 
values of V which must be considered. The general trends 
are discussed in §6. 
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5.3.2 Correlation Coefficient Analysis 

We considered the Spearman partial rank correlation co- 
efficients between the f our relevant so urce characteristics 
P,D,z and a. Following Macklin (1982), we calculated the 
partial rank correlation coefficients with four variables, e.g., 



rpD.za = 



rpD,z — r Pa ,r Da 



1/2 > 



(10) 



for the correlation between P and D independent of z and 
a. Here the three-variable partial correlation coefficient is 



tpd,z = 



rpp — rp z rp z 
(1 - r Dz f (1 - r Pz ) 2 



(11) 



with rpD being the Spearman correlation coefficient between 
two variables P and D. 

The significance level associated with the 4- variable cor- 
relation is 

EpDtZa = (jW~5) 1/2 ln ( l+r PD , za , ^ 2) 



1 - r PD , za J ' 

where, N samp is the size of the sample considered. 

The relevant Spearman partial rank correlation coeffi- 
cients involving [P, D, z, a] for the data and for the models 
for those cases for which the 2-D KS tests were done are 
given in Table 7. The four-variable correlation coefficients 
(rpD, Z a,rp z ,Da, etc) were computed by combining the ob- 
served data or the model "simulated" data for all the rele- 
vant surveys: 3C, 6C and 7C III. We do so in order to dilute 
the tight [P — z] correlation in a single flux-limited com- 
plete sample (BRW), and to detect the correlations which 
exist between other characteristics. The 2-variable correla- 
tion, rpD, was always negative; however, in the 4-variable 
case, with the effects of z and a removed, rpo, Z a showed 
a small positive correlation. We also examined the correla- 
tion coefficients of the data and model simulations in each 
survey, 3C, 6C, or 7C, separately. 



6 DISCUSSION 

We have performed quantitative tests of three detailed mod- 
els for RG evolution. This is the first attempt to perform 
such statistical tests involving 4 radio source observables 
over 3 complete radio surveys. During our multi-dimensional 
Monte Carlo simulation procedure we found that it is very 
difficult to get acceptable simultaneous fits to the radio prop- 
erties P, D, z and a for all three redshift complete subsam- 
ples of the 3C, 6C and 7C radio catalogs. This is true using 
either the default parameters suggested by each of these 
three leading models, or when considering extensive varia- 
tions upon them involving changing one or more of the pa- 
rameters to plausible different values. Usually the P and z 
fits were correlated, due to flux limiting arguments discussed 
before. The fits to the 6C survey were generally better com- 
pared to those for 3C and 7C; this is because of the smaller 
number of sources in the 6C catalog and the nature of the 
KS test. Our weighting of the "total 1-D KS probability" by 
the square root of the number of sources helps to compen- 
sate for this. It was most difficult to get acceptable fits to 
the faintest sources cataloged in 7C. 

When varying the model parameters from their default 
values the greatest improvement came from steepening the 



power law index for the initial beam power d i stribu tion to 
x = 3 from x = 2.6 used bv iBlundell et alJ <1999h . This 
change improved the KDA and MK model performances 
greatly (Tables 3 and 5). The KS statistics for the BRW 
models were never wonderful, especially the D fits (Table 
4); nonetheless, varying the maximum age assumed for the 
sources from 500 Myr to 150 Myr for the KDA and MK 
models and to 250 Myr for the BRW model also produced 
better fits. 

We found the following trends for the ratios of number 
of sources detected in the 6C and 7C simulations and the 
number in the actual catalogs, as compared to 3C (Ratioec 
and Ratio7c respectively). For the KDA and BRW models, 
the detection number ratio was more consistent for 6C than 
for 7C simulations; i.e., Ratioec was closer to 1.0 (which it 
should equal ideally) than was Ratio7c . For the MK model, 
the detection number ratios for 6C and 7C (which were in 
the range 0.7— 1.2) were equally consistent. Though we cal- 
culated the detection number ratios, we do not display them 
nor did we formally consider them in comparing the models. 
These ratios can be made closer to 1 by varying the redshift 
birth function or the RLF (Eq. 1), and so are not good tests 
of the models per se. 

From the 2-D KS test results we find that the [P — z], 
[P — D] and [z — D] planes can be reasonably fitted by the 
"varied" models, particularly those for KDA and MK. Most 
of those probabilities are > 0.2 for the KDA and most exceed 
0.05 for MK. All of the P's of the "varied" BRW model 
not involving a are higher than those of the default BRW 
model. Improvements are also seen for all of the non-a MK 
P's using the "varied" model. This is the case for only 7 of 9 
Vs of the KDA "varied" model. These models cannot fit any 
plane involving a, with all the a-related 2-D KS probabilities 
^ 0.01 for every model. These 2-D results provide support 
for the hypothesis that the "varied" models based on 1-D 
KS tests are indeed better fits. 

By comparing the values of 2-D KS probabilities in the 
models of Table 6, we conclude that KDA model is the best 
(having the highest number of "P's close to 1) in fitting the 
observational data, very closely followed by MK, and finally 
BRW. 

From the 4-variable correlation coefficient results (Ta- 
ble 7) we see that the KDA model is able to match the sur- 
vey data correlations very closely (at least for P, D, z). The 
matches to the data correlations are less good for the BRW 
and MK models. The parameter variation cases which were 
the best fits (i.e., gave highest combined P[p y D, z ,a}) when 
judged with respect to 1-D and 2-D KS tests, are not neces- 
sarily the better cases according to the correlation analyses. 
The KDA default performs better than the KDA varied (1-D 
KS best fit) case. For BRW and MK models, the default and 
the varied cases perform comparably (i.e., sometimes the de- 
fault verison is a better match to the data correlations and 
sometimes the varied fit is better). 

Considering the signs of the four- variable coefficients, 
the MK model predicts [P — a] anti-correlation and [D — z] 
correlation which are trends opposite to the survey data and 
to the other models. The sign of the [D — a] correlation of 
the surveys is only predicted by MK, while the other models 
produce an anti-correlation; however, given the very poor a 
distribution for the MK model this advantage is meaning- 
less. 
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Figure 2. The [P — D — z — a] planes for the observational samples 3CRR, 6CE and 7CRS. The symbols classify the sources into redshift 
bins as follows; Plus: ^ z < 0.5, Triangle: 0.5 ^ z < 1.0, Cross: 1.0 ^ z < 1.5, Square: 1.5 ^ z. The P — z correlations, arising from the 
flux limits, are clear in the first row. The D — z plane shows the decreasing trend of average size as redshift increases. 



From the correlation coefficient analyses we conclude 
that the KDA model fits the data most closely, followed 
by BRW, and finally MK. Similar trends are also seen if we 
examine the coefficients obtained by considering each survey 
separately. 

We plotted slices through the [P — D — z — a] volume (P 
vs z, P vs D, D vs z, and a vs z) for each of the simulated 
surveys, and examined their consistency by comparing them 
with the overall trends in the [P — D — z — a] planes of 
the actual data. The actual data is shown in Fig. 2. The 
simulated data are shown in Figs 3, 4 and 5 for the KDA, 
BRW and MK models, respectively. The plotted simulations 
are for one of the best (in a statistical sense) parameter sets 
for each model. Plots for other good parameter values appear 
similar, while those for worse parameters (according to our 
KS summary statistic) look less like the data. Sources are 
detected out to similar values of redshift, power and size in 



the 3C simulations as in the data. The KDA and MK models 
show very similar trends in P, D, and z. The unique features 
of the BRW model results are discussed below. 

Unsurprisingly, the values of P and z exist in a cluster 
in the P — z plane, above a lower curve determined by the 
flux limit of the survey. All of our simulated surveys of all 
models miss many of the low z - low P sources seen in the 
data. Very high z sources (z > 2.5) are underproduced in 
all the 7C simulations, and a similar, but less pronounced, 
trend is also present for 6C. All the 3C simulations present 
a greater scatter in P for high P - high z sources (z > 1) 
when compared to the data. A few powerful, high z sources 
are detected in the 3C simulations at z > 2.0 which are not 
present in the data. The scatter in P is naturally less in the 
6C survey because of the upper (as well as lower) flux limit. 

Examining the P — D planes of the simulations, we 
find that the KDA and MK models overproduce small and 
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Figure 3. The [P — D — z — a] planes for the 3C, 6C and 7C simulations of the KDA Model. The initial ensemble is formed using 
x = 3.0, TMaxAge = 150 Myr; the power evolution is with parameter changes po = Po (Default) /2 = 3.6 X 10~ 22 kg m~ 3 , p m = 2.12, for 
a case with initial source population size = 4861474. The symbols are as in Fig. 2. 



large high power sources in 3C, and underproduce the large 
weaker sources. The underproduction of low z sources is 
manifested in the P — D planes of the 6C and 7C simu- 
lations as the absence of less powerful sources (due to the 
P — z correlation). 

There is a strong P — D evolution seen in the BRW 
model, which is most pronounced in the 3C and 6C sim- 
ulations. The 3C simulation overproduces powerful smaller 
sources and misses several large ones. The 6C and 7C simula- 
tions underproduce less powerful, smaller sources. Again, the 
KDA and MK models show too weak P — D anti-correlations 
than does the data (at least for 3C), whereas the BRW model 
shows too strong an anti-correlation. 

The 6C and 7C simulations show a paucity of low z and 
high z sources in the D — z planes of all the models. The 
KDA and MK models overproduce very small and very large 
3C sources at all redshifts. The BRW simulation presents a 



stronger anti-correlation of linear size with redshift, specially 
for 3C, where there are no large sources at intermediate 
redshifts. 

The D— z evolution (decrease of D as z increases) occurs 
due to imposing survey flux limits. This is a ramification of 
the "youth -redshift degeneracy" discussed in §5.1. The high 
redshift sources show a very steep decline of their luminosi- 
ties with age (seen from the P — D tracks in Fig. 1) and fall 
below the survey flux limits at young ages, as their radiat- 
ing particles undergo severe inverse Compton losses off the 
CMB and adiabatic expansion losses as they are transported 
from the high pressure hotspot to the lobes. Thus, we can 
only detect these high z sources at an early age when they 
are still above the limiting survey flux. These younger high z 
sources are naturally smaller and yield the weak "linear size 
evolution" (seen in the D — z plane). Both the KDA and 
MK simulations do not show this effect as clearly as does 
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Figure 4. The [P — D — z — a] planes for the 3C, 6C and 7C simulations of the BRW Model. The initial ensemble is formed using 
x = 3.0, TMaxAge = 250 Myr; the power evolution is with parameter change an = 7.5 kpc, for a case with initial source population size 
= 3355926. The symbols are as in Fig. 2. The upward arrow in the 3C panels implies that one data point exists outside the plotted range 
of the figure. 



the actual data. On the other hand, the BRW simulations 
show stronger D — z anti-correlations than do the data. 

There are several observational features (including 
trends in the [P — D — z — a] planes of the data samples) 
that cannot be explained by any models considered so far. 
The [P — D] diagram for the 3CRR data show a clear anti- 
correlation with large scatter. Another interesting feature 
is the clump of sources in the 6CE [P — D] diagram nea r 
D ~ 100 kpc, Pisi ~ 27.5 W Hz" 1 sr" 1 iNeeser et alll99Eft . 
Neither of these is reproduced in the models. The KDA and 
MK model simulations predict too many very large D > 1 
Mpc and powerful sources (more in 3C, some in 7C), which 
ar e not present in the data. T his feature has been discussed 
in iKaiser fc AlexanderlJl999D . 

The BRW V(D) were very low for many cases (espe- 
cially for 3C), and the BRW [P — D] diagrams for all 3 sim- 



ulated surveys showed too strong a [P — D] anti-correlation. 
This arises because the BRW model simulations produce too 
many small but powerful sources. A possible explanation of 
this problem could be synchrotron self-absorption of the ra- 
diation emitted by such small powerful sources, which is not 
included in the model. Thus some small sources should fall 
below the survey flux limit at a frequency of 151 MHz. In- 
cluding this effect could improve the relative performance of 
the BRW model. 

An important point to remember is that all three mod- 
els considered here are incomplete in the sense that they 
do not incorporate enough physics to predict the complete 
physical conditions prevailing in FR II radio sources. 

Consideration of additional factors may be necessary in 
these models. First, the environmental density (p) could vary 
with redshift and it must eventually deviate from its power 
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Figure 5. The [P — D — z — a] planes for the 3C, 6C and 7C simulations of the MK Model. The initial ensemble is formed using x = 3.0, 
TMaxAge = 150 Myr; the power evolution is with parameter change ^ ma x(hs) = 3 X 10 s , for a case with initial source population size = 
4861474. The symbols are as in Fig. 2. 



law behavior with distance. The beam power (Qo) distri- 
bution might vary with redshift and the maximum lifetime 
of AGN activity (TMaxAge) could vary with redshift and jet 
power. Also, the birth function of radio sources with redshift 
(RLF), could have a greater variation with luminosity. 



7 CONCLUSIONS AND FUTURE WORK 

We have compared the leading models of radio lobe power 
evolution for FR II RGs, namely the KDA, BRW and MK 
models, using a simulated radio survey prescription (follow- 
ing BRW). Each of the dozens of simulated radio surveys 
we computed required the generation and analysis of a few 
10 6 to > 10 7 radio sources and hence substantial amounts of 
computing power. The total number of Monte Carlo simu- 
lations done exceeded 250 and over a billion individual RGs 



were evolved; this was necessary to narrow down the set of 
parameters for each model to the "best fit" ranges described 
in the present work. One-dimensional KS tests were done to 
narrow down the parameters of the different models to lo- 
cate more desirable ones. These preferred parameter sets of 
the models were then compared with the data by using 2-D 
KS tests, and correlation coefficient analyses. 

Hydrodynamical modeling of classical dou bl e ra- 
dio sources (e.g . - [llooda. Mangalam fc Wiital Il994l ; 
ICarvalho fc O'Deal I2002T) shows that the pressure in 
the nearly self-similarly growing lobes falls wit h time while 
the ho tspot pressure does not vary much. The lKaiser et all 
( 1997) model examined here assumed that the head pres- 
sure falls with time (and is proportional to that of the 
cocoon), so this is a weakness of that model. iBlundell et alJ 
( 1999) adopted a constant hotspot pressure (implying more 
adiabatic losses for particles in the hotspot of older sources) 
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while considering the adiabatic expansion of particles out of 
the hotspots to the lobes. They showed a rough qualitative 
agreement between thei r simulated and real 3C an d 7C data 
in the [P — D — z] space. iManolakou fc K irk (2002) modified 
the BRW picture by proposing an acceleration mechanism 
occurring throughout the head region; they obtained [P—D] 
tracks in somewhat better accord with 3CRR data, but did 
not consider [P — D — z — a] distributions. 

Our much more extensive simulations and statistical 
analyses, based on KS tests and correlation coefficients, pro- 
vides a quantitative way to directly compare these three 
models. We note that despite the hundreds of simulations 
we computed which did employ substantial variations on 
the default parameters for each RG model (only a portion of 
which are displayed in this paper) we could not completely 
cover the entire plausible parameter space. We also note 
that other figures of merit could have been devised to distin- 
guish between the goodness of fits of the data to the various 
simulation results, since no really suitable multi-parameter 
statistic is available for samples of this size. Keeping these 
caveats in mind, we believe both that we have covered the 
vast majority of the sensible parameter ranges and that our 
choice of combined KS probabilities is a good way to com- 
pare different simulations. In this spirit, we now present our 
conclusions comparing how the models performed in differ- 
ent aspects of consistency between the simulations and data. 

Our key result is somewhat disappointing. Despite in- 
vestigating a wide range of parameters we find that no exist- 
ing model gives excellent fits to all the data simultaneously. 
However, from the statistical test results, the KDA model 
appears to give better fits than do the BRW or MK models. 

Explicitly judging from the 1-D KS test results, the MK 
model frequently produces acceptable statistics for P, z and 
D. The KDA simulations also often give adequate statistics. 
The BRW simulations do not give as good statistics as do the 
MK and KDA models. After incorporating the 10 kpc linear 
size cut-off the statistics for some BRW models improve, but 
are still not as good as those given by the other two models. 

According ot the 2-D KS test results, the KDA model 
fits the data most closely, then comes MK, and finally BRW. 
From both the 1-D and 2-D KS test results, planes in the 
[P — z — D] space can be reasonably fitted by some parameter 
variation of the models (fits determined by higher values of 
the KS probabilities), but it is difficult to get acceptable 
a fits. Both the KS tests dictate that the "varied" model 
parameters of all the models (Tables 6 and 7) are better fits 
to the data as compared to the default parameter values. 

However, in terms of reproducing the correlations be- 
tween the source properties (Spearman partial rank corre- 
lation coefficient), the default models perform better than 
(KDA) or comparable to (BRW, MK) the "varied" models. 
The KDA model correlations match the survey data corre- 
lations most closely, followed by BRW and then MK. 

Our ana lyses used the redshi ft birth function of radio 
sources from Iwillott et al.l teOOllTs radio luminosi ty func- 
tion. We conclude that, using IWillott et al.l (1200 ill 's RLF, 
the KDA and MK models perform better than BRW in fit- 
ting the 3CRR, 6CE and 7CRS survey data when compared 
with respect to KS-based statistical tests, and the KDA 
model provides the best fits to the correlation coefficients. 

This is the first in a series of papers which aim at com- 
paring the performances of radio source evolution models. 



Our goal is to develop one which is a good fit to all the ob- 
servational data. We are performing similar tests on a mod- 
ified model we are developing, whose results will be pre- 
sented elsewhere. This new model incorporates conical jet 
expansion for a fraction of a radio source's lifetime within 
the BRW and MK models. This allows us to incorporate a 
variable hots pot size in the models, wh ich is supported by 
observations dJevakumar fc Sa ikia 2000|). Here the hotspot 
pressure varies as a function of the linear size or the source 
age, which is a more practical possibility for hotspots of RGs 
evolving over 100s of Myr. 

In the future, we plan to extend this work by allow- 
ing redshift variations in the environmental density profile 
(in particular, we will allow for variations of po, ao and (3 
with cosmic epoch). We also will consider jet propagation 
through ambient media which change from power law den- 
sity decline to constan t densities (which change with z) at 
scales around 100 kpc. iBarai et alJ (120041) gives preliminary 
work on the implications of the volumes attained by radio 
sources considering cosmological evolution of the ambient 
gas density. 

Our final aim is to estimate the volume fraction of the 
relevant universe occupied by radio lobes, and hence to test 
the robustness of the exciting, but preliminary, conclusion 
that expanding radio galaxies play a significant role in the 
cosmological history of the universe. 
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Table 3. KDA Model: 1-D KS Statistics for Selected Parameter Variations 1 
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Table 3 — Continued 



X Model P[P,D,z,a] 

T MaxAg e h Ensemble Size Survey V(P) V(D) V{z) V(a) P[p, 2 D, z , a ] 
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4.89e-17 


2.16 


150 


4861474 


6C 


0.879 


0.244 


0.308 


1.29e-09 


2.72 






70 


0.0271 


0.626 


0.0544 


3.39e-08 




3.0 


po = pi c 


3C 


0.413 


0.0502 


0.766 


7.42e-21 


2.33 


150 


p = 2.12 


60 


0.879 


0.113 


0.293 


1.29e-09 


2.60 




4861474 


70 


0.105 


0.232 


0.0544 


3.24e-07 





a These results do not exclude sources with total linear size D(t) < 1 kpc. 
h T M axAge in units of Myr. 

c Parameter values set equal to those given in the first KDA model jKaiser et alJll997f) . 

d Parameters defining external environment density profile are set to those of the BRW model, 
namely 13 = 1.6, a = 10 kpc, p = 1.67 x 10" 23 kg m -3 . 

C P1 = PO (Default)/2 = 3.6 X 10~ 22 kg m - 3 . 
f P 2 = 5 X po (Default) = 36 x 10~ 22 kg m" 3 . 
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Table 4. BRW Model: 1-D KS Statistics for Selected Parameter Variations' 1 



X 


Model 












7~^[P,D ,z ,a] 


rp b 

J- AI ax Age 


Ensemble Size 


Survey 


V(P) 


P(D) 


Viz) 


V(a) 


W\ POD is- rJ\ 


Z.O 


Default 


oL 


4.yue-n 


l.z4e-lU 


1.4oe-Uo 


n noon 
U.Uozz 


U.U0D4 


500 


2930490 


6C 


1.93e-05 


0.0240 


0.00120 


3.66e-10 


0.0728 






70 


1.94e-09 


1.20e-04 


1.96e-07 


0.0115 




2.0 


Deiault 


OL 


l.lze-lU 


i .ole-lo 


i at^ nn 

1.4 1 e-uy 


n c\a ceo 
U.U40U 


u.tyt 


250 


1466378 


6C 


0.00415 


9.92e-04 


0.189 


3.66e-10 


0.203 






7C 


5.03e-08 


0.00101 


9.96e-09 


0.0198 




o.U 


JJeiault 




o.oye-Uo 


t.yye-ty 


i kc;^ n/i 


C OPC 

o.t>ye-Uo 


O AOQQ 

U.Uzoo 


500 


2930490 


6C 


5.06e-04 


9.50e-05 


0.0162 


3.66e-10 


0.0241 






7C 


4.17e-07 


3.70e-05 


8.07e-07 


0.0194 




o.U 


JJeiault 




8.o4e-Ui> 


i.yue-iy 


n nm i o 
U.UUllz 


1 A K« C\A 

i.4oe-U4 


n one 
U.2U8 


250 


1571349 


6C 


0.00789 


2.13e-04 


0.277 


3.66e-10 


0.229 






7C 


1.66e-06 


1.22e-04 


2.18e-05 


0.0104 




o . U 


P — 1 


oL 


f\ C\C\A SO 

U.UU4oz 


n nnt^/io 
U.UUo4y 


U.U44o 


1.44e-U / 


n c^nv 
U.oU / 


250 


1571349 


6C 


0.0893 


0.113 


0.214 


3.03e-17 


0.772 






70 


2.75e-05 


0.306 


1.56e-04 


1.13e-06 




o . U 


do — ZU Kpc 


oL 


l.o ( e-U4 


A OQn 1 O 

4.Zo6-lZ 


n nni 
U.UUloO 


i n^o c\a 
l.Uoe-U4 


n Q^ii 
U.ODl 


250 


1571349 


6C 


0.0283 


0.00415 


0.214 


8.69e-14 


0.566 






7C 


1.01e-06 


0.333 


7.58e-06 


2.26e-05 




o.U 


Po - Pi 


■J f 1 

6\j 


1 7o« nt; 
1. / ze-Uo 


l.loe-lo 


U.UU28 1 


l.Uze-U I 


n oon 

u.zzy 


250 


1571349 


6C 


0.0432 


0.00756 


0.214 


1.59e-12 


0.292 






7C 


6.62e-08 


0.0957 


2.08e-06 


6.91e-04 




o.U 


t bf = 1UU yr 


oL 


o.ole-Uo 


b.U4e-ty 


n nm i o 
U.UUllz 


A (1(1^ OCC 

4.yye-uo 


n 1 oc 
U.loo 


250 


1571349 


6C 


0.00405 


4.62e-04 


0.277 


3.66e-10 


0.192 






7C 


4.25e-06 


0.0103 


7.36e-06 


0.0237 




3.0 


p = 2.001 


3C 


0.264 


5.03e-06 


0.680 


2.87e-04 


1.39 


250 


3355926 


6C 


0.0473 


0.511 


0.141 


3.66e-10 


1.72 






7C 


9.07e-05 


0.00538 


8.04e-06 


0.00332 




3.0 


t;, s = 10 3 Yr 


3C 


0.240 


1.10e-07 


0.474 


0.00487 


1.42 


250 


3355926 


60 


0.0481 


0.509 


0.310 


3.66e-10 


1.89 






7C 


7.53e-04 


0.251 


1.46e-04 


0.00332 




3.0 


ao = 7.5 kpc 


3C 


0.330 


2.90e-09 


0.668 


9.34e-07 


1.63 


250 


3355926 


6C 


0.268 


0.00764 


0.583 


1.29e-09 


1.73 






7C 


1.55e-05 


0.148 


5.04e-05 


7.21e-04 





a These results are without any total linear size cut-off. 
b T Ma xA 3 e in units of Myr. 

C A11 other parameters are as given in the BRW model l)Blundell et alJll99flh . 

V = 2 x po (Default) = 3.34 x 10~ 23 kg m- 3 . 
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Table 5. MK Model: 1-D KS Statistics for Selected Parameter Variations' 1 



x Model V[p,D,z, a ] 

Ensemble Size Survey V(P) ~P(D) V{z) V(a) P[p,2£>,z,a] 



2.6 


Default 


3C 


3.77e-12 


2.68e-05 


3.99e-07 





0.270 


500 


4397469 


6C 


0.0283 


0.0685 


0.308 


1.80e-24 


0.324 






7C 


5 38e-04 


0.0119 


0.0199 


2 67e-18 




2.6 


Default 


3C 


4.05e-14 


0.00375 


4.97e-06 





0.670 


150 


3888492 


6C 


0.0149 


0.254 


0.585 


1.80e-24 


0.961 






7C 


9 77e-09 

•J • lit; \J\J 


0.413 


1 18e-05 


2 1fie-15 




3.0 


Default 


3C 


0.00819 


0.0356 


0.229 





1.23 


500 


11236430 


6C 


0.0419 


0.353 


0.580 


1.54e-18 


1.83 






7C 


6 56e-04 


0.541 


00708 


2 lfie-15 




3.0 


KDA Env. d 


3C 


0.126 


4.01e-04 


0.353 





0.842 


500 


/3, oo,po 


6C 


0.212 


0.135 


0.134 


3.14e-23 


0.935 




1 4659499 


7C 


001 08 


0.0149 


068Q 


6 Q0p-1 7 




3.0 


P = 1.6 


3C 


0.267 


0.288 


0.351 





1.73 


500 


14659422 


6C 


0.150 


0.311 


0.169 


7.86e-21 


2.38 






7C 


0593 


0.647 


0.176 


9 1 6p-1 5 




3.0 


do = 7.5 kpc 


3C 


0.282 


0.0741 


0.106 





1.76 


500 


14659422 


6C 


0.117 


0.674 


0.438 


1.43e-17 


2.61 






7C 


n os^q 

u.uooy 


809 


0889 


9 (S^p 17 




3.0 


Po = pi ° 


3C 


0.400 


0.0913 


0.130 





1.63 


500 


14659422 


6C 


0.0770 


0.205 


0.335 


1.54e-18 


2.30 






7C 


0839 


0.756 


0889 


9 65p-1 7 




3.0 


^fmin(hs) 7 


3C 


0.170 


0.0961 


0.147 


1.12e-44 


1.44 


500 


14659422 


6C 


0.361 


0.624 


0.376 


3.54e-16 


2.00 






7C 


0.0210 


0.265 


0.0794 


7.04e-21 




3.0 


"fmax(hs) = 3 X 10 8 


3C 


0.290 


0.0692 


0.171 





1.53 


500 


14659422 


6C 


0.404 


0.0540 


0.321 


8.73e-23 


1.98 






7C 


0685 


0.519 


0.132 


9 1 6p-1 5 




o.U 


— 9 nni 

p — Z.UUl 




n HG1 7 


o.Doe-uo 






1 IK 

l.oO 


500 


9772948 


6C 


0.819 


0.00394 


0.183 


2.60e-20 


1.48 






7C 


0.00750 


0.00908 


0.0970 


1.94e-08 




3.0 


e=1.5 


3C 


0.00323 


0.211 


0.702 





1.96 


500 


11236430 


6C 


0.0468 


0.258 


0.618 


3.19e-18 


2.78 






7C 


4.48e-04 


0.704 


0.00779 


2.76e-18 




3.0 


Default 


3C 


0.0277 


0.420 


0.633 





1.95 


150 


4861474 


6C 


0.310 


0.420 


0.182 


3.14e-23 


2.92 






7C 


6.75e-04 


0.267 


0.0247 


2.16e-15 




3.0 


po = pi e 


3C 


0.161 


0.0186 


0.766 





1.97 


150 


4861474 


6C 


0.879 


0.0135 


0.0790 


1.80e-24 


2.39 






7C 


0.00807 


0.638 


0.0183 


1.91e-13 
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Table 5 — Continued 



X Model V[P,D,z,at] 

T MaxAs e h Ensemble Size Survey V(P) P(D) V(z) P{a) V{ P , 2 D,z, a ] 



3.0 


13 = 1.6 


3C 


0.0473 


0.0406 


0.808 





2.37 


150 


4861474 


6C 


0.351 


0.780 


0.394 


2.71e-20 


3.40 






7C 


0.0238 


0.942 


0.0795 


1.58e-15 




3.0 


7m«i(ft») = 3 x 10 8 


3C 


0.0678 


0.0361 


0.668 





2.47 


150 


4861474 


6C 


0.740 


0.809 


0.199 


1.80e-24 


3.58 






7C 


0.00815 


0.917 


0.0267 


1.54e-15 





a These results do not exclude sources with total linear size D(t) < 1 kpc. 
b TMaxA 3 e in units of Myr. 

c P arameter values set equal to as given in the preferred MK model (case B, iManolakou fc Kirkl 

d Parameters of external medium density profile are set to those of the default KDA model, namely, 
P = 1.9, a = 2 kpc, p = 7.2 x 10~ 22 kg m -3 . 

°pi = pa (Default)/!- 5 = 1.133 x 10" 23 kg m" 3 . 
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Table 6. 2-D KS Test Results for the Three Models 



Model 
Parameters 


Survey 


P(P - 2) 

' V / 


V(P - D) 


2-D KS Probability, P(K-S) 
V(z-D) P{P-a) 


Viz - a) 


V(D - a) 


J\iJA 




i .uoe-UD 


4.z i e-uy 


3.99e-07 


4.79e-09 


7 cc c no 


1 A Kcs HQ 


Default a 


6C 


0.816 


0.0164 


0.00741 


3 38e-05 


4.71e-05 


2.17e-04 




7C 


0.0108 


0.0124 


0.00876 


0.0135 


0.00763 


0.00192 






U.ool 


U.UZOo 


0.129 


5.70e-15 


4.0DC- 14 


( .ooe-zo 


Varied b 


6C 


0.445 


0.370 


0.244 


9 88e-04 


9.15e-04 


0.00458 




7C 


0.00832 


0.251 


0.226 


8.73e-05 


3.21e-05 


5.51e-07 


D rCW 




1.4/ e-Uo 


q on,, i n 
o.zUe-lU 


1.06e-08 


1.22e-08 


z.Uoe-Uo 


o.yUe-14 


Default a 


6C 


3.08e-04 


0.00116 


n nnQ/i/i 


z.oic-uy 


2.81e-08 


5.87e-05 




7C 


5.89e-08 


7.40e-07 


3.15e-06 


2.04e-04 


3.09e-05 


4.30e-06 


BRW 


3C 


0.432 


1.34e-07 


5.70e-07 


3.68e-05 


3.36e-05 


2.71e-17 


Varied c 


6C 


0.654 


0.00918 


0.0205 


1.19e-05 


2.92e-05 


1.81e-05 




7C 


2.30e-04 


5.89e-04 


1.60e-04 


8.75e-04 


1.29e-04 


5.06e-05 


MK 


3C 


5.05e-10 


3.20e-13 


3.40e-09 


4.09e-35 


1.55e-35 


3.28e-31 


Default a 


6C 


0.0843 


0.0491 


0.221 


2.24e-16 


3.88e-17 


3.27e-13 




7C 


1.98e-04 


4.04e-04 


4.81e-03 


4.72e-10 


2.36e-10 


9.57e-09 


MK 


3C 


0.0431 


0.0117 


0.0701 


2.24e-35 


2.73e-35 


8.38e-34 


Varied d 


6C 


0.177 


0.510 


0.244 


2.75e-17 


9.82e-18 


1.50e-13 




7C 


0.0247 


0.0688 


0.0885 


3.72e-10 


1.87e-10 


1.05e-08 



Simulations with the respective model parameters as used by the authors. Initial population generated 
using x = 2.6, T M axA ge = 500 Myr. 

b KDA model simulation using initial population with x = 3.0, TMaxAge = 150 Myr. The power evolution 
is with parameter changes po = po (Defauit)/2 = 3.6 x 10 -22 kg m -3 and p = 2.12, other parameters set to 
their default values, for a case with initial source population size = 4861474 (last parameter variation entry 
of Table 3). 

C BRW model simulation using initial population with x = 3.0, TMaxAge = 250 Myr. The power evolution 
is with parameter change ao = 7.5 kpc, other parameters set to their default values, for a case with initial 
source population size = 3355926 (last parameter variation entry of Table 4). 

d MK model simulation using initial population with x = 3.0, TMaxAge = 150 Myr. The power evolution is 
with parameter change r y ma x(hs) = 3x 10 s , other parameters set to their default values, for a case with initial 
source population size = 4861474 (last parameter variation entry of Table 5). 
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Table 7. 4-variable Spearman Partial Rank Correlation Analysis a 



Data Model (combining all surveys a ) 

KDA BRW MK 



Coeff. 


All a 


Default 


Varied b 


Default 


Varied b 


Default 


Varied b 


rpD,za ° 


0.0303 


0.0528 


0.198 


0.102 


0.0944 


0.358 


0.196 


V d 

^PD,za 


0.478 


0.844 


3.20 


1.63 


1.51 


5.96 


3.16 


Tp z ,Da 


0.716 


0.668 


0.648 


0.415 


0.576 


0.569 


0.495 


~Spz,Da 


14.2 


12.9 


12.3 


7.04 


10.4 


10.3 


8.63 


TDz,Pa 


-0.268 


-0.274 


-0.206 


-0.106 


-0.234 


0.303 


0.433 


^Dz,Pa 


-4.33 


-4.48 


-3.33 


-1.70 


-3.78 


4.97 


7.37 


rp a ,Dz 


0.147 


0.0456 


0.318 


0.329 


0.428 


-0.167 


-0.103 


Sp a ,Dz 


2.33 


0.729 


5.25 


5.45 


7.27 


-2.68 


-1.65 


TDa,Pz 


0.472 


-0.0287 


-0.640 


-0.890 


-0.881 


0.922 


0.888 




8.08 


-0.459 


-12.1 


-22.7 


-22.0 


25.5 


22.5 


r za ,PD 


-0.0234 


0.0970 


-0.0935 


-0.0237 


-0.226 


-0.465 


-0.569 


^za,PD 


-0.369 


1.55 


-1.50 


-0.379 


-3.65 


-8.01 


-10.3 



a The four observables P, D, z and a for the 3C, 6C and 7C III surveys (whether real 
or simulated), combined together in a single sample. 

b The particular parameters used are the same as those in Table 6 for each of the KDA, 
BRW and MK varied models. 

c Spearman partial rank correlation coefficient between two variables P and D, when 
the other two variables z and a are kept fixed. 

Significance level associated with the correlation between P and D, independent of 
z and a. 



